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We investigate the phase behaviour of 2D mixtures of bi-functional and three-functional patchy particles and 
3D mixtures of bi-functional and tetra-functional patchy particles by means of Monte Carlo simulations and 
Wertheim theory. We start by computing the critical points of the pure systems and then we investigate how 
the critical parameters change upon lowering the temperature. We extend the Successive Umbrella Sampling 
method to mixtures to make it possible to extract information about the phase behaviour of the system at a 
fixed temperature for the whole range of densities and compositions of interest. 



I. INTRODUCTION 

The role of the particle valence — defined as the ability 
to form only a limited number of bonds with neighbour- 
ing particles — in controlling the phase behaviour of col- 
loidal systems has been emphasized in numerous recent 
studie^i^. Limited valence has emerged as the key ele- 
ment in the formation of equilibrium colloidal networks, 
commonly named gels. In this state of matter, rigidity 
is enforced by the long-life of the inter-particle bonds, 
which, at low temperatures T, is longer than the exper- 
imental observation time^^. Lowering the valence guar- 
antees that these low T gel states do not clash with phase 
separation, which is always present when the thermal en- 
ergy becomes significantly lower than the characteristic 
bond strength. In contrast with the out-of-equilibrium 
colloidal gels formed via spinodal decomposition followed 
by a kinetic arrest induced by the strong depletion attrac- 
tion^, limited valence gels are equilibrium states. It has 
indeed been shown that lowering the valence shifts the 
phase separation boundaries to low densities^ ^, opening 
a wide region of particle concentrations where stable gels 
can fornP. 

Binary mixtures of limited valence particles enrich con- 
siderably more the spectrum of possibilities offered by 
limited valence. First, the average valence of the system 
can take non-integer values, which has been exploited to 
investigate the approach to the limit of valence two. In 
this limiting case, phase separation is completely sup- 
pressed, since particles aggregate in long chains that in- 
teract only through excluded volume interactions. Sec- 
ond, mixtures can be exploited to tune the selectivity of 
the network to different species, stabilizing mixed or in- 



terpenetrating gels^^ ^^, the equilibrium equivalent of the 
recently reported out-of-equilibrium bigels^^. 

In limited valence systems, gas-liquid phase separation 
arises from a subtle competition between the number of 
bonds that can form in the two coexisting phases (the en- 
ergy term) and the entropy. The latter accounts for the 
degeneracy of the bonding patterns, which is different in 
the two phases. The gas phase is usually formed by di- 
luted clusters, while the liquid phase is characterized by 
a percolating network of bonds. This competition is cap- 
tured by the thermodynamic perturbation theory (TPT) 
developed by WertheimP^Hi^ to model the behaviour of 
associating fiuids, the atomic and molecular analogues of 
limited- valence colloids. Wertheim theory is a powerful 
tool for investigating the ph ase b ehaviour of pure fiuids 
as well as of binary mixtures^^EHl^ 

Recent theoretical studies of binary mixtures with dif- 
ferent compositions have revealed a subtle interplay be- 
tween the entropy of bonding and the entropy of mixing, 
with a marked effect on the phase diagram of the mix- 
ture. Interestingly, Wertheim theory predicts that for 
binary mixtures of bi- and three- functional particles, the 
gas-liquid critical density pc and the critical tempera- 
ture Tc decrease as the average valence decreases {i.e., 
on increasing the fraction of bi-functional particles/^, in 
agreement with existing numerical results^. By contrast, 
for a binary mixtures of bi- and tetra-functional particles, 
Wertheim theory predicts a qualitative difference in the 
behaviour of the critical parameters. Indeed, binary mix- 
tures of such species are expected to initially undergo an 
increase of pc as the average valence decreases. A further 
increase in the number of bi-functional particles inverts 
this trend and pc decreases and tends to a constant value 



as the fraction of bi- functional particles approaches one. 
Finally, mixtures of two and five functional particles are 
predicted to have critical densities which increase mono- 
tonically as the average valence is lowered^^. 

The numerical evaluation of phase coexistence in bi- 
nary mixtures is not an easy task. Indeed, the presence 
of a second species adds a new axis to the phase diagram, 
which is now a three-dimensional volume defined by T 
and the density of the two components (T — pA— Pb) or 
by other combinations of pA and pB as T — p — x, where 
P = Pa -^ pB is the total number density and x = pb/p 
is the composition of the mixtur d^^ * ^ . The two coex- 
isting phases are characterized by different values of x. 
Even focusing on a specific x value (along the so-called 
dilution line), the determination of the phase diagram 
requires the evaluation of the shadow lines^^. In this 
manuscript we introduce a new and powerful computa- 
tional method to investigate the phase behaviour of mix- 
tures in the whole three-dimensional volume, by extend- 
ing the successive umbrella sampling method^, which 
was shown to be very effective in the evaluation of the 
phase behaviour of single component systems^^^^. By 
applying this new methodology it is possible to compute 
the entire density of states for the binary mixture, i.e., 
the information required to estimate phase coexistence at 
all compositions. As a test case of scientific relevance we 
compute the phase diagram of mixtures of bi-functional 
and three-functional patchy particles in two-dimensions 
(21^) and of bi-functional and tetra-functional patchy 
particles in three dimensions {3D). We then compare 
numerical results with theoretical predictions based on 
Wertheim theory, confirming the predicted growth of the 
critical density in the 2 — 4 mixture as the fraction of 
bi-functional particles is increased. The analysis of the 
calculated phase behaviour shows that the pc growth on 
increasing the number of bi-functional particles results 
from a progressive transformation of the transition from 
condensation to demixing. At the same time, increasing 
the number of particles with two patches does reduce the 
region in density where the instability takes place, con- 
firming the general trend that a reduction of the average 
valence increases the density region where a stable gel 
can form. 



II. METHODS 



connecting the centre of particle i with the patch a on 
its surface. The function / can then be written as 



f{^^.^J) 



1 if 



^ii • Vf > COS I' 



otherwise, 



Vj > cos On 



for any a, 

for any /3, 



A. Model 



(1) 

where ^max controls the width of the patches. All the 
patches share the same shape, i.e., 6 and 6>max are fixed 
and do not depend on particle species. 

The only difference between particles of different 
species is in the number of patches patterning their sur- 
faces. In 2D we study binary mixtures of particles deco- 
rated by either 3 (species A) or 2 (species B) patches^^. 
For particles of species A, the patches are symmetrically 
placed on the equator, for particles of species B they are 
located on the poles. The KF parameters are S = 0.03a 
and cos^max = 0.894. 

In 3D we study a binary mixture with the two species 
having either 4 (species A) or 2 (species B) patches. 
Patches are located on the particle surface in a tetra- 
hedral fashion^^ for species- A particles, on the poles for 
species-5 particles. The KF parameters are S = O.llQa 
and cos^max = 0.92. 

Both sets of parameters fulfill the geometrical single- 
bond-per-patch condition sin6>max < 1/(2(1 + ^)), Pre- 
venting patches from being involved in more than one 
bond^. 



B. Computational methods 

1. Pure systems 

To compute the location of the (pseudo-)critical points 
of the pure systems we rely on the Bruce- Wilding (BW) 
mixed-field scaling method^^. This technique provides 
an expression for the order parameter M which can be 
used to fit the probability distribution P{M) of non- 
symmetric fiuids to the Ising one, in order to estimate 
the pseudo-critical parameters of the finite system. On 
top of that, the BW approach provides scaling expres- 
sions which can be used to obtain the values of the crit- 
ical parameters in the thermodynamic limit. With this 
procedure, the deviation from the average value of the 
order parameter M at criticality can be written as^ 



Each particle is modelled as a hard sphere (in 3D) or 
a hard disk (in 2D) of diameter <j, decorated with a fixed 
number of interacting patches on the surface. The patch- 
patch interaction between particles i and j is described 
by a Kern-Frenkel (KF) potential, i.e., it is a square- 
well potential of range S and depth e, modulated by a 
function f{fti^flj) which depends solely on the particle 
orientations Cti and Ctj. Let f^j be the normalized vector 
joining the centres of particles i and j and vf the versor 



AM = M - Mc oc p^su, 



(2) 



where p = N/V is the number density, N is the number 
of particles, V is the volume of the system, u = U/V is 
the energy density and 5 is a non-universal {i.e., model- 
dependent) factor. 

In order to calculate the joint probability distribution 
p{N^U)^ to be compared with the probability distribu- 
tion of the Ising order parameter, we rely on Grand 



Canonical Monte Carlo (GMMC) simulations, i.e., sim- 
ulations at fixed temperature T, number of particles TV 
and chemical potential ja. We employ successive umbrella 
sampling (SUS)^^ to overcome the high free-energy barri- 
ers between the two phases. With this method, the region 
to be explored is partitioned in overlapping windows of 
AN particles. Each window is then sampled with GCMC 
simulations with appropriate boundary conditions'^, pro- 
viding a speed-up proportional to the number of windows 
explored in parallel. 

To properly locate the pseudo-critical point we make 
use of Eq. ^ to project p{N, U) and to obtain p(AM). 
In turn, we extract the critical parameters by matching 
this distributio n to the (2D or 3D) Ising order parame- 
ter distribution'^ '' by means of histogram reweighting 
techniques^. 



2. Mixtures 

Now we introduce this extended SUS method for the 
case of a generic mixture. Let S be the number of species 
and [Nf^^'^.,Nf^^^] be the range of number of particles 
of species i of interest. Applying the SUS method con- 
sists in partitioning the [Ar™,Arfa^] x [N^''' , N^''''] x 
... X [A/'^^", 7V^^^] space into overlapping windows of size 
ni X n2 X . . . X n^. Without any loss of generality we can 
fix the width of the overlap to be 5w and then the total 
number of windows is n^ = ni=i \{Nf^^^ — Nf^^^)/{ni — 
5w)'\^ where [•] stands for the ceiling function. Each 
window is identified by a iS-dimensional index w. All 
the windows are then explored through special GCMC 
simulations, i.e., simulations performed at fixed T, V, 
{/i^}, with the additional constraints that the number 



(a) 



of particles of species i, A^i, has to lie within the range 
[N^ .,N^ + rii] for each i G [1,5']. The main simulation 
output is the histogram counting how many times a state 
with A^i , A'2 , . . . , A^s- has been visited in a given window 
w^ namely p^(A'i, A'2, • • • , A's'). The total free-energy 
density profile p(A'i, A'2, . . . , A's-) is then computed by 
using the overlapping portions of the windows to join to- 
gether all the p^ . This operation is done through a least- 
squares method. Let p^(A^i, A^2, • • • , A^s) be the partial, 
already joined part of the total p(A'i, A'2, . . . , A'^) and 
p^(A'i, A'2, . . . , A^s-) be the histogram of the next win- 
dow to be attached. Then, in order to extend p^ to the 
{A"^} values stored in p^ one needs to multiply the latter 
by the factor h^ given by 



(3) 
where {O^} is the set of A"!, A'2, ... , Ns values which are 
in the overlapping region between pP and p^ . 

There is no unique way of performing this operation, 
since the process of joining different windows can start 
from any window and follow any pattern. For the present 
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FIG. 1. Description of the scheme employed to reconstruct 
the complete density of states p{Na, Nb) for a binary mixture. 
(a) Raw output from the nine 3x3 windows with the lowest 
number of particles. The arrows show the order with which 
the windows are joined together, (b) The same data has been 
used to compute the final p{NajNb)'- there are no visible 
boundaries between the windows. The brighter the color, the 
higher the value of p. The curves have been smoothed out in 
order to increase readability. 



study the scheme delivering the best results is the follow- 
ing. Since the system under study is a binary mixture, 
we compute the joint probability distribution p(A^A, A"^), 
i.e., we keep track of how many times the system has A"^ 
particles of type A and A"^ particles of type B. We start 
with the (0, 0) window and then begin to attach windows 
along the first species' direction, so that the second win- 
dow is (1,0), the third is (2,0) and so on. Once the last 
window in the row has been joined, we start attaching an- 
other row of windows by joining the (0, 1) window. This 
procedure is schematically shown in Figure [T] 

At low temperatures, i.e., when the numerical noise 
increases, using Eq. (|3| on the simulation data may re- 
sult in histograms which cannot be reliably reweighted 
at all the required chemical potentials. This happens be- 
cause the main contributions to the factors b'^ are the 



largest p^{Na^ Nb) values and, since reweighting to very 
different chemical potentials moves the signal to the less 
precisely-attached windows, the quality of the resulting 
histogram deteriorates. To overcome this difficulty we 
first reweight each p^ and then join them together by 
means of Eq. (|3| . 

The simulation output p{Na-,Nb)^ computed at fixed 
T, jiA and /i^, can be evaluated at different chemical 
potentials /i^, fi'^ by histogram reweighting, i.e., 

(4) 
where /3 = l/(/cT), here k is the Boltzmann constant. 
Therefore, p{Na-, Nb) encodes all the information on the 
system in the whole investigated Na^ Nb plane at fixed 



where c is a fitting parameter which depends on temper- 
ature. Note that, unlike to what we do in pure systems, 
we do not store any information on the energy of the 
system, and hence we do not perform any temperature 
reweighting. 

In the rest of the article, the superscript ^^ (^^) is used 
to refer to quantities associated to the 2D {3D) model. 

The simulation box sizes L^^ = 16.9 and L^^ = 10 
are kept fixed throughout this work. We do not perform 
any finite-size scaling study and hence we compute only 
pseudo-critical parameters. For the sake of brevity, in the 
following we use critical instead of pseudo- critical when 
referring to these quantities. 



C. Theory 



In order to obtain information on the phase behaviour 
of the mixture we employ the following criterion: we 
reweight p{Na^Nb) at a certain jua and then we tune 
jiiB until p{Na^Nb) is double-peaked, with the area be- 
low the two peaks being equal. If no such fiB value can 
be found, then we are out of the coexisting region for 
the chosen (T^jha) values. If the equal area condition 
is fulfilled, the total free-energy density profile can be 
split up as a sum of two contributions Pi{Na,Nb) and 
P2{Na^Nb)^ one for each phase. In the systems studied 
here, this is done by making a cut in 7V^, Nb plane at 
fixed Na = N^^ where N^ is the position of the fitted 

minimum of the p{Na) = '^n^=iP{Na^ Nb) curve. The 
number of particles of species i in the phase j, (N^) is 
then computed by taking an average over the appropriate 
particle number distribution, i.e., 



(Ni) 



EnI=oEn!=Ori,P^iNA,NB) 

Armax Armax 



(5) 



The quantities can, in turn, be used to compute com- 
positions x^^^ = (A^B V((^A^ + ^B^)) and densities 



The pressure P can be computed by considering that 
p(0, 0) = e~^^^ is the grand-canonical partition func- 
tion, and hence^ 



kT 



log(p(0,0)). 



(6) 



Similarly to what we do for pure systems, we com- 
pute the (pseudo-) critical points of mixtures by compar- 
ing the p(M), obtained by projecting the p{Na^Nb)^ to 
the Isin g orde r-parameter distribution of the right dimen- 
sionality^^^, with the only difference being the choice of 
the order parameter^, defined as 



AM (X Pa ^ cpB, 



(7) 



We also investigate the patchy colloidal mixture theo- 
retically by means of Wertheim's ffist order perturbation 
theory. A detailed description of the original theory can 
be found in Refs.^^^^^. Here we briefly quote the results 
and set the not ation fo r Wertheim's theory extended to 
binary mixturea^^l^^^. The Helmholtz free energy per 
particle of the mixture is: 



Jh = Fh/N = fref + h 



(8) 



where N = Na-\-Nb is the total number of particles, fref 
is the free energy per particle of the reference fluid of hard 
spheres (HSs) in 3D or hard disks (HDs) in 2D, and /& is 
the bonding free energy per particle. As usual we write 
fref as the sum of ideal-gas and excess terms: fref = 
fid + fex- The ideal-gas free energy is given (exactly) by 

f3fid{v.x^'^)=\nr]-l^ ^ x^'hn{x^'^Vi), (9) 

i=A,B 

where Vi is the (irrelevant) thermal volume, x^*^ = Ni/N 
is the molar fraction of species i = {A, 5}, and r] = 
Va + Vb is the total packing fraction {r] = VsP, with p 
the total number density and Vg = 7r/6<j^ the volume of 
a HS in 3D or Vs = n/Aa'^ the area of a HD in 2D). The 
excess part accounts for the excluded volume interactions 
between the monomers. Both species have the same size 
and hence we can approximate the excess part by the 
well-know Carnahan-Starling equation of state for hard 
spheres in the 3D mixture: 



PfeAv) = ^^4§ (3D) , 



(10) 



(l-r?)2 

and use the HendersoiP^ equation of state for hard disks 
in the 2D mixtures: 



PfeAv) = -lHi-v) + ^ '' 



8(1-??) 



(2D) . (11) 



The bonding free energy is approximated by Wertheim's 
thermodynamic first-order perturbation theory 



/3/fc = {M) InX. 



Nu 



(12) 



where X„ is the probability that one site is not bonded 
and 



(M) = a;('4)M('4)+a;(^)M(^), 



(13) 



is the average number of patches per particle in the 
mixture (M^^^ and M^^^ are the number of patches of 
species A and B respectively). The probability of finding 
an unbonded patch is related to the total density, molar 
fractions and absolute temperature through the law of 
mass action 



X„ = 1-7?X2a„(M) 



(14) 



The bond between two patches is characterized by A^^. 
Using the Kern-Frenkel potential, we find: 



An 



1 



^(r) [exp(/3e) - 1] dr. 



(15) 



where ^(r) is the radial distribution function of the refer- 
ence fluid of HS or HD, and the integral is calculated over 
the volume (area) of a bond, v^. If v^ is small enough, 
we can approximate the radial distribution function by 
its contact value, gdv). Under this assumption the Eq. 
( 15 ) simplifies to 



An = -gM [exp(/3e) - 1] . (16) 



The contact value of the radial distribution function is 

1-^/2 



9civ) 
for hard spheres and 

9c{v) 



(l-r?3) 



(3D) 



V 



1-7? 16(l-7?)2 



(2D) 



(17) 



(18) 



for hard disks. The bonding volume (area) is related 
to the depth and range of the potential. For the three 
dimensional mixtures we find 

Vb = 7r{l-cos0maxf{{cr^Sf-cr)/3 = 0.00269cr^ (3D), 

(19) 
and for the two dimensional case 



Vb 



92 

^max 



([a + Sf - a'^) = 0.00418(7^ (2D). (20) 



Finally, we obtain the equilibrium properties of the mix- 
ture by minimising (at a fixed pressure, composition 
and temperature) the Gibbs free energy per particle 
gcix^ p, P, T) = P/p + fn with respect to the total den- 
sity. We use a standard Newton-Raphson method to 
minimise go- Coexisting points are located by a stan- 
dard common-tangent construction on gcix) at constant 
temperature and pressure. Critical points are computed 
by determining those states which satisfy the spinodal 



condition, fwfxx — {fxv)'^ = 0- In addition, stability re- 
quires the vanishing of the third-order derivative in the 
direction of largest growth: 



Jvvv "^Jx 



'^JxVV \ r j Jvvv \ r I — U, 

\ Jvv J \ Jvv J 

(21) 

where subscripts denote partial derivatives, i.e., j^v is 
the second partial derivative of fn with respect to the 
reduced volume per particle v = I/t^ and the composition 
X (molar fraction of bi-functional particles) at constant 
temperature. 



III. RESULTS 
A. Critical points 

In both investigated binary mixtures, particles of 
species B are bi-functional and therefore do not exhibit 
any gas-liquid phase separation^ ^^. On the other hand, 
particles of species A have a higher valence and exhibit a 
regular gas-liquid phase separation at low densities and 
temperatures^^. 

We start by computing the pseudo-critical parame- 
ters of the pure systems, i.e., of the tetra- functional 
model in 3D and of the three- functional model in 2D. 
We then move down in temperature and compute the 
phase boundaries on the whole Na^ Nb plane at fixed 
temperature. Figure |2] shows the best fits to the 2D and 
3D Ising order parameter distribution for all the binary 
systems studied in this work. Data are less scattered in 
the 2D case but, despite the noise, the resulting fits are 
reliable at all temperatures. 

From the p{Na,Nb) reweighted at criticality we can 
extract all the critical parameters of the mixtures, such as 
critical composition Xc, critical density pc, critical pres- 
sure Pc and average valence (M). Figure [3] shows the 
critical parameters of all the computed critical points ac- 
cording to Monte Carlo simulation. In a mixture of two- 
and three-functional patchy particles in 2D, see panel (a), 
the critical density goes down as the critical temperature 
decreases. Moreover, pc is a monotonically increasing 
function of (M), as shown in panel (b). By contrast, the 
mixtures of two- and tetra-functional patchy particles in 
3D display a pc drop as the average valence is increased, 
see panel (d), as predicted by Wertheim's theor}^^^ xhis 
growth is reflected in the behaviour of pc versus Tc, panel 
(c), which is monotonically decreasing in the investigated 
range of temperatures. In both systems the dependence 
on Tc of the critical pressure Pc (not shown) follows that 
of the critical density in this range of temperatures: Pc 
increases with Tc in 2D and decreases in 3D. 

The results from Wertheim's theory are depicted in 
Figure [4J The theory predicts the same behaviour as 
Monte Carlo simulation for all the critical parameters 
and both types of mixtures considered here. The agree- 
ment is quantitative only for the critical temperature and 
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FIG. 2. Best fits to the Ising order parameter distribution 
for mixtures of two- and three-functional patchy particles in 
2D (a) and two- and tetra-functional patchy particles in 3D 
(b). Open symbols: Monte Carlo simulation at different re- 
duced temperatures T* = kT /e. Solid black line: Ising order 
parameter distribution. 



the average valence at the critical point. Wertheim's first 
order perturbation theory does not include the formation 
of closed loops of patchy particles. As a result, the crit- 
ical density is underestimated. This discrepancy may 
also arise from the different nature of the critical phe- 
nomenon, which is mean-field in Wertheim's theory and 
Ising in simulations. 

In 2D (bi- and three-functional particles) and in the 
limit of zero temperature (note that using Wertheim's 
theory we can investigate the whole range of tempera- 
tures) the critical density vanishes, panel (a), and the 
valence at the critical point tends asymptotically to two, 
panel (b). Mixtures of bi- and tetra-functional particles 
in 3D behave differently, see panels (c) and (d). The the- 
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FIG. 3. Critical parameters for (a)-(b) a mixture of hard disks 
with two and three patches and for (c)-(d) a mixture of hard 
spheres with two and four patches, according to Monte Carlo 
simulation. Panels (a) and (c): reduced critical temperature 
as a function of the critical density. Panels (b) and (d): crit- 
ical density as a function of the average number of patches 
per particle at the critical point. 



ory predicts a reentrant behaviour for the critical density 
occurring for kT/e < 0.12, i.e., for temperatures which 
are lower than those studied here by Monte Carlo sim- 
ulations. At very low temperatures the critical density 
tends asymptotically to a value different from zero, and 
the average valence at the critical point is always higher 
than two. The critical pressure (not shown) also shows 
a non-monotonic behaviour upon lowering T. It first in- 
creases near the critical temperature of the pure tetra- 
functional fluid, and then decreases. For mixtures of bl- 
and five-functional particles (not shown) in 2D and in 
3D the critical density increases monotonically as T is 
lower ed^^. 

To conclude this section on critical points, we com- 
pute the critical composition Xc, i.e., the fraction of bi- 
functional particles (species B) at criticality, as a func- 
tion of Tc. This quantity, shown in Figure [5J is monoton- 
ically decreasing for both types of mixtures, as predicted 
by the theory. Therefore, as the temperature is lowered, 
the composition of the fluid at criticality tends to favour 
more and more the bi- functional particles. At very low 
temperatures Xc tends asymptotically to one in mixtures 
of bi- and three- functional particles (2D), and to ^ 0.9 
in mixtures of two- and tetra-functional particles (3D). 



B. Coexistence region 

We proceed to analyse the phase boundaries. For the 
mixtures considered here, upon lowering the tempera- 
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FIG. 4. Critical parameters for (a)-(b) a mixture of hard disks 
with two and three patches and (c)-(d) for a mixture of hard 
spheres with two and four patches, according to Wertheim's 
theory. Panels (a) and (c): reduced critical temperature as a 
function of the critical density. Panels (b) and (d): critical 
density as a function of the average number of patches per 
particle at the critical point. The empty circles indicate the 
position of the critical point for a pure fluid of three- (tetra-) 
functional particles in 2D (3D). The region explored via MC 
simulations is shaded in grey. 



ture the instability region is first encountered at the crit- 
ical temperature of the pure system, namely Tc . If we 
project the three-dimensional phase diagram onto the p^, 
pB plane, then the phase-coexisting region is a point ly- 
ing on the pB = ^ axis at T = Tc . If T decreases, the 
instability region expands towards larger values of pB^ in 
line with the results for the critical parameters reported 
in the previous section. This T-dependence is shown in 
Figure [6) which displays the phase boundaries and the 
computed critical points projected onto the pA^ Pb plane 
for both types of mixtures. The temperature dependence 
in the explored T range is qualitatively different in the 
two cases. 

The overall density of the low-density phase in the 
two dimensional 2 — 3 mixture, see panel (a), does not 
change much upon changing T. The high-density phase 
branch, on the other hand, displays a rather strong T- 
dependence: as T is lowered, it extends to larger values 
of p2 and smaller values of p^. Therefore, at all the inves- 
tigated temperatures the phase transition has the char- 
acter of an ordinary gas-liquid phase separation, with 
a rather broad separation in densities between the two 
phases. 

Three dimensional 2 — 4 mixtures, however, behave in 
a different manner, see panel (b). Even though the effect 
of T on the liquid branch is similar to that observed in 
the 2D system, the phase originating from the gas-phase 
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FIG. 5. Reduced critical temperature, kTc/e^ as a function of 
composition (i.e. fraction of bi-functional particles) at the 
critical point, Xc, for mixtures of bi- and three- functional 
patchy particles in 2D and bi- and tetra-functional patchy 
particles in 3D. (a) Results from Monte Carlo simulation, (b) 
Results according to Wertheim's theory. The region investi- 
gated with MC simulation is shaded in grey. 



of the pure system does not remain confined at low den- 
sities but it moves at larger p2 and p4. Indeed, at the 
lowest temperature, the difference in overall density be- 
tween the two coexisting phases is remarkably smaller 
than in the 2D case. In this case the phase transition 
has a stronger demixing component which tends to segre- 
gate the bi-functional particles in the low-density phase. 
This may be ascribed to a change in the balance between 
the entropy of mixing and the entropy of bonding as the 
number of patches n on particles of species A changes. 
When two particles of species A form a bond, the result 
is a two-particle cluster with 2(n — 1) sites available for 
bonding. By contrast, a cluster of two particles of species 
A and B has n available sites. Therefore, the difference 
between the number of available bonding sites on the two 
types of clusters is n — 2 sites. That is, as n increases, 
forming a cluster of dissimilar species is less favourable 
entropically than having a cluster made up of particles 
of the same species. The gain in the entropy of mixing, 
on the other hand, remains the same. As a result, the 
tendency for phase separation increases with n. Ref.^ 
contains a more detailed analysis on this topic. 
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FIG. 6. Projection of the phase diagram of the (a) 2D binary 
mixture in the p2, p3 plane, and (b) 3D binary mixture in the 
y92, pA plane according to MC simulation. Symbols: densities 
of the coexisting phases for different reduced temperatures 
T* = kT/e. The symbols connected by the dashed line mark 
the positions of critical points. The orange filled circle marks 
the position of the critical point associated to the pure system 
of three (four) patches in 2D (3D), occurring at Tc — 0.136 
(Tc — 0.157). Note that tie lines connecting two coexisting 
points (not shown) are not vertical lines in this plane. 



The same results according to Wertheim's theory are 
presented in Figure [7[ As expected, the theory underes- 
timates the coexisting densities both in 2D and 3D. As 
shown in the previous section, the 2D critical line tends 
asymptotically to p2 = and pa = 0, which makes it pos- 
sible to find coexisting phases at arbitrarily low densities. 
By contrast, the 3D critical line ends at finite densities 
p2 and /)4. In addition, in the 2D mixture the theory 
captures the overall behaviour of the system. In partic- 
ular, the decrease of p2 and pa at criticality, observed in 
simulations, is well reproduced. In the 3D case, on the 
other hand, there is a qualitative difference between nu- 
merical and theoretical results. As the system is cooled 
down, there is an increase of the critical p2 and p4 as 
computed in simulations, while the theory predicts a de- 
crease of the critical /)4. This difference is associated to a 
different slope of tie lines, i.e. of the lines connecting the 
two coexisting phases. Tie lines for the lowest-T systems 
are shown in Figure ^(simulations data) and in Figure^ 



FIG. 7. Projection of the phase diagram of the (a) 2D binary 
mixture in the y92, ps plane, and (b) 3D binary mixture in the 
P2, pA plane (b) according to Wertheim's theory for different 
reduced temperatures T* = kT/e. The black dashed line is 
the line of critical points of the mixture. Tie lines are not 
vertical lines in this representation. 



(theoretical results). 

If tie lines have positive slopes, the coexisting points 
follow the same trend as constant-composition lines 
which, in this representation, are straight lines with zero 
intercept. This means that high-density phases always 
have a higher density of both species than the coexist- 
ing low-density phases. In the systems studied here, this 
is true for the 2D numerical and theoretical results and 
for the theoretical 3D mixture. By contrast, the numeri- 
cal 3D system exhibits tie lines with negative slope: the 
gas-like phases always have more bi-functional particles 
than their respective liquid-like phases. In this regard, 
the theory fails to capture the demixing nature of the 
phase separation occurring in the 3D system as observed 
in simulations. On passing, we note that Wertheim the- 
ory predicts tie lines with negative slopes in mixtures of 
particles with two and six patches (now shown). 

A different representation can be constructed by pro- 
jecting the three dimensional phase diagram onto the p, x 
plane, i.e., total density against composition (molar frac- 
tion of bi-functional particles). This is shown in Figure 10 
(Monte Carlo simulation) and in Figure 11 (Wertheim's 



theory). This representation makes it clear that the in- 
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FIG. 8. Tie lines for the lowest-temperature system in (a) 2D 
and (b) 3D (dashed lines). Orange diamonds pinpoint criti- 
cal points. The shadow-cloud construction is highlighted in 
panel(a). The dash-dotted line is the dilution line at fixed 
composition x — 0.181. It crosses the coexistence region 
in two points which are part of the cloud curve (filled red 
squares). Their coexisting partners are part of the shadow 
cloud (open red squares). 



crease of the critical density in the 3D systems is due to 
the contribution of the low-density phase, which moves 
to larger densities and compositions as the temperature 
is lowered in the vicinity of the critical temperature of 
the pure four-patches fluid. In the 2D system, on the 
other hand, the density range of the low-density phase 
does not change strongly with T, while the composition 
increases steadily. This results in the shift of the critical 
density to smaller values. 

This particular projection of the phase diagram hides 
some of the differences we previously noted between the- 
ory and simulations. Indeed, tie lines in this representa- 
tion always have negative slope. 



C. Cloud and shadow curves 

A common bi-dimensional representation of the phase 
diagram of binary mixtures is done by making a cut at 
a fixed composition in the three-dimensional pA-,pB-,T 
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FIG. 9. Theoretical results for tie lines at the lowest- 
temperature system in (a) 2D and (b) 3D (dashed lines). 
Open circles mark the positions of critical points. 



phase diagram. The intersection between the cutting 
plane and a fixed temperature plane is usually called di- 
lution line^^, because by following it the density of the 
system can be varied without changing the overall com- 
position. The points at which the dilution line intersect 
the phase boundary are the so-called cloud points. Each 
of these points coexists with a infinitesimal amount of 
the other phase which has, in general, a different com- 
position. These points are called shadow points. An 
exemplification of this procedure is given in Fig. [SJ^a). 
The shadow and cloud curves are then constructed by 
plotting the sets of cloud and shadow points on the T, 
p plane. Note that this cloud-shadow construction is a 
projection and therefore, while the cloud points have the 
same composition, the composition of the shadow points 
varies, i.e., the tie lines are out-of-plane. The spinodal 
lies inside the cloud curve and the two lines touch at the 
critical point, which is always located at the intersection 
of the shadow and cloud curves^^^^. 



Figure [12] shows the cloud-shadow construction for the 
2D system at a fixed composition x = 0.181. The or- 
ange diamond marks the position of the critical point at 
kT /e = 0.12, which happens to have the same composi- 
tion. In line with the results of Section |IIIB[ the curves 
resemble the usual bell-shaped gas-liquid coexisting re- 
gions observed in pure systems^ ^^ ^^. Both cloud and 
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FIG. 10. Projection of the phase diagram of the (a) 2D bi- 
nary mixture and the (b) 2D binary mixture in the p, x plane 
according to MC simulation. Symbols: densities of the coex- 
isting phases for different reduced temperatures T* = kT/e. 
The symbols connected by the dashed line mark the positions 
of critical points. The orange point marks the position of the 
critical point of the associated pure system with three (four) 
patches in 2D (3D). Note that the tie lines are not vertical 
lines. 



shadow curves have a low-density (high-density) phase 
which is in coexistence with a high-density (low-density) 
phase. In addition, all the coexisting branches have a 
monotonic T-dependence. Indeed, both the low-density 
cloud branch and the high-density shadow branch have 
nearly T-independent densities, whereas the densities of 
the low-density shadow branch and of the high-density 
cloud branch are monotonically decreasing and increas- 
ing, respectively, as in gas-liquid phase transitions of 
pure, simple fluids. All the cloud-shadow constructions 
computed for the 2D mixtures share the same qualitative 
features with the one shown here, regardless of x. The- 
oretical results, shown in Fig. 12 ^b), are qualitatively in 
line with simulation data. 



The 4 — 2 mixture, as previously noted, behaves in a 
different way. Figure [13] shows the cloud and shadow 
curves for a mixture with a fixed composition x = 0.09. 
The whole cloud curve has a monotonic T-dependence 
and its behaviour is qualitatively similar to the 2D case. 
The shadow curve, on the other hand, exhibits several 



FIG. 11. Projection of the phase diagram of the (a) 3D binary 
mixture and the (b) 2D binary mixture in the />, x plane 
for different reduced temperatures T* = kT/e according to 
Wertheim's theory. The black-dashed line is the line of critical 
points. 



remarkable features. First of all, the density of the low- 
density branch is rather high, being half-way between 
the two branches of the cloud curve and very close to the 
critical density. Moreover, its temperature dependence 
is non-monotonic, first decreasing and increasing upon 
lowering T. The high-density branch seems to be non- 
monotonic as well and, at a slightly lower temperature 
than the lowest investigated T, a crossing between the 
cloud and the shadow high-density branches is expected. 
We ascribe the lacking of this feature in the theoretical 
curves, shown in Fig. 13 ^b), to the demixing character of 
the transition. 



IV. CONCLUSIONS 

We have extended the successive umbrella sampling 
method^^ to binary mixtures in oder to simulate the full 
bulk phase diagram. The method consists in dividing the 
simulation cell into a set of overlapping windows of vari- 
able size, which is set by the number of particles of each 
species allowed in a particular window. Each window is 
then sampled by means of grand canonical Monte Carlo 
simulations, counting the number of times that different 
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FIG. 12. (a) Cloud (black circles) and shadow (green squares) 
curves for the 2D binary system at fixed x — 0.181. Each 
point on the gas branch of the cloud curve (filled circles) is in 
coexistence with the point on the liquid branch of the shadow 
curve at the same temperature (filled squares) , and vice versa 
(empty symbols) . Lines with arrows point out the connection 
between cloud and shadow points, (b) Theoretical results for 
the same system at the same value of x. The orange diamonds 
signal the position of the critical point of the mixture having 
the same composition. 



microstates, characterized by the number of particles of 
each species, appears in each window. This gives the 
probability distribution of microstates, and thus the de- 
sired thermodynamic properties of the system. We have 
tested the validity of the new method by computing the 
phase behaviour of patchy colloidal mixtures. We have 
compared the results with the theoretical predictions of 
Wertheim's first order perturbation theory. 

We have analysed two types of colloidal mixtures: bl- 
and three- functional patchy particles in two dimensions, 
and bi- and tetra-functional patchy particles in three 
dimensions. In the first case, bi- and three- functional 
patchy colloids in 2D, simulation and theory are in excel- 
lent agreement. The agreement is quantitative for those 
variables that do not involve the density, such as the 
temperature or composition. The theory underestimates 
the density (a well know problem of Wertheim's first or- 
der perturbation theory) as it neglects the formation of 
closed loops. As a consequence the agreement between 
the predicted densities and those found in the simulations 
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FIG. 13. (a) Cloud (black circles) and shadow (green squares) 
curves for the 3D binary system at fixed x — 0.09. The symbol 
filling has the same meaning as in Fig. |12[ The orange dia- 
mond signals the position of the critical point of the mixture 
having the same composition, extrapolated from the location 
of the nearby computed critical points, (b) Theoretical results 
for the same system at the same value of x. 



is only qualitative. The phase behaviour of these mix- 
tures is similar to that of b i- and three-functional parti- 
cles in three dimension^Ell; ^j^e critical density vanishes 
as the critical temperature approaches zero, suppressing 
condensation and yielding an increasingly large region of 
phase space where empty liquids are stable. Thus, di- 
mensionality does not change the topology of the phase 
diagram of these mixtures. This was to be expected 
since the energy and entropy of bonding dominate the 
behaviour of such mixtures and this is determined by the 
functionality of the particles, rather than by the spatial 
dimension. 

Three dimensional mixtures of bi- and tetra-functional 
colloids, however, exhibit a different behaviour. The re- 
sults of Wertheim's theory predict that as the critical 
temperature decreases the critical density first increases, 
then decreases and tends asymptotically to a value, which 
is different from zero. The simulation confirms the ini- 
tial increase of the critical density but the non-monotonic 
behaviour could not be confirmed by simulations as it oc- 
curs (according to the theory) at very low temperatures, 
a region not accessible by the current simulation tech- 
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niques. Despite the qualitative description of the mix- 
ture's critical behaviour, the theory fails to describe ad- 
equately the shadow and cloud curves for these mixture. 
The origin of this discrepancy may be related to the for- 
mation of closed loops of particles, neglected by the the- 
ory, which increase in number as the functionality of the 
particles increases. 

In summary we have developed and tested a new simu- 
lation scheme to investigate the phase diagram of binary 
mixtures of patchy colloidal particles. The new method 
can be applied to a large variety of problems, such as 
the surface and confinement properties of patchy col- 
loidal mixtures, or the study of the bulk and percolation 
properties of the recently predicted bicontinuous gels or 
bigel j^^tti^i ^ i^ mixtures of patchy particles these bigels 
may be equilibrium structures, when they occur in the 
empty liquid regime^ or dynamically arrested structures 
when they occur inside the liquid-vapour or the liquid- 
liquid binodals, as in ordinary binary mixtures^^. A de- 
tailed investigation of the connectivity and other physical 
properties of these structures, in and out of equilibrium, 
is bound to reveal novel features with potential applica- 
tions. 
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